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NUMERICAL INTEGRATION OF ORBITS IN MULTIREVOLUTION STEPS* 


by 

C. E. Velez 
Goddard Space Flight Center 


INTRODUCTION 

Many artificial satellite orbit determination problems, such as lifetime 
studies or computing long-range predictions, require the numerical integration 
of the equations of motion over many periods of revolution. An algorithm de- 
signed to save computing time in solving such problems has been proposed in a 
paper by Cohen and Hubbard (1960). Essentially, the method involves combining 
the numerical integration with a procedure which "steps" the calculations ahead 
in multirevolution increments. This "stepping" procedure, both in formulation 
and usage, is similar to the Adams predictor -corrector process in numerical 
integration. Briefly, a cycle in the algorithm consists of first extrapolating or 
"predicting" the orbital elements n-revolutions ahead. Then, starting with these 
extrapolated values, the equations of motion are integrated over one revolution. 
Finally, the extrapolated values are improved by using a corrector formula along 
with the information obtained from the single -revolution integration. 

The multirevolution "predictor" or extrapolation formula was introduced in 
a paper by Mace and Thomas (1960). In this paper it was demonstrated that if 
successive revolutions are similar, so that the successive descending nodes are 
close together, one may extrapolate without loss of accuracy as many revolutions 
as will make the change in the orbital elements from node to extrapolated node 
about the same as the change in step-by-step integration over one revolution. 
Supplementing this work, Cohen and Hubbard developed a corrector formula for 
this n-revolution extrapolation. In addition, it was also shown that this corrector 
formula reduces to the Adams formula in the limiting case, i.e. , the coefficients 
of the corrector were shown to be polynomials in 1/n with the constant terms 
being the Adams coefficients. 

In this paper, the similarity in derivation between the multirevolution for- 
mulas and the Adams predictor corrector formulas will be explored and used 
to develop recursive formulas for the coefficients. In addition, the multirevo- 
lution algorithm will be detailed, and examples demonstrating its usage will be 
given. 

*This report is a republication of NASA X-5 14-67-341, January 1967. 
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FORMULATION 


Let f be a function defined on an interval [x Q , x m ] with equal spacing h (i.e. , 
x. -x._ 1 = h, i = 1 , 2, • • • m), and assume f is continuously differentiable. 
Consider the difference operators V, A, E a , D, and 1 defined as 


and 


Vf(x) = 

f(x) - f(x-h) 

(backward difference), 

Af(x) = 

f (x + h) - f (x) 

(forward difference^ 

E a f(x) = 

f ( x + ah ) 

(shifting operator), 

Df(x) = 

f' (X) 

(differentiating operator), 

1 f(x) = 

f(x) 

(identity) , 

xe [ x o> X J- 

Also, if L is any one of these operators, define L n by 


L n f(x) = L[L n_1 f ( x ) . 


Some basic relations between these operators are* 


A = 


V 


(l) 


E = 1 + A , (2) 

and 

E = e hD . (3) 


We begin by deriving the Adams integration formulas. The predictor can 
be derived by expressing the forward difference operator in terms of backward 

•Equality interpreted in the sense of operators. For verification of these formulas and justification for the 
algebraic manipulation of these operators, see Hildebrand (1956), pp. 128-134. 
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differences of the derivative. We have from (2) and (3) that 


hD = - log (1 -V) 


( 4 ) 


Hence, using (1), 


V V 1 

A “ l-v - 1 - V - log ( 1 - V) 1x0 




v.i= 0 


r D. 


(5) 


where the coefficients a. can be determined in the following manner, 

l 

Using the expressions 


and 


1 - x 


x 


CO 


i=0 


(6) 


- log (1 - x) 


00 

V" 1 X 1 
x + 1 ’ 


i=0 


(7) 


we have 


V 1 

1 - V - log ( 1 - V) 


CO 

L 

i=0 


Z ■ r ~ T ~ 1 


00 

= 2_V*. 

i=0 


i = 0 
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where the a. are given recursively by 


a o “ 1 


and 


a. - 1 


■ L a a i-j 

j=l 


( 8 ) 


i = 1, 2, 3, 




1 


where 

= 

Hence, using (5) and (8), our formula emerges in the form 
Af(x) = f(x +h) - f(x) 


( 9 ) 


h^l+ 2 ^ V 2 + • • f ' ( x ) 


( 10 ) 


which is the familiar Adams -Bashf or th formula used to obtain a "predicted" 
value of f (x + h) when solving first-order initial -value problems. The corre- 
sponding corrector formula can be derived by expressing the backward differ- 
ence operator in terms of backward differences of the derivative. 

As above, using (4), 


V 

V " -log(l-V) hD 



where the coefficients a.* are given recursively by 

*o* = 1 

and 


CL* ~ 


1 





( 11 ) 


( 12 ) 


where the / 3 - ’s are given by (9). 
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Then, using (11) and (12), 


Vf(x) - f(x) - f(x-h) 

“ , (13) 

which is the Adams-Moulton formula used to obtain a corrected value of f(x) in 
the numerical integration process. 

From (5) and (11), 

£l + V + V 2 + * • • J + a.j* V + • • • J = a Q + a.j V + • • • , 


so equating coefficients of like powers, we have the relationship 

k 

a* 

1 

i=0 



(14) 


We now wish to derive the multirevolution predictor and corrector formulaSc 
This is done precisely as above, except that differences of the derivative are re- 
placed by differences at a large interval of differences at a small interval. To 
this end, we define the differences A and V by 

’ n n J 


and 


A n f ( x ) " f(x +nh) - f(x) 
V „ f ( x ) = f ( x ) ~ f ( x ~ nh ) • 


where n is same fixed integer, i.e., A^ and are the forward and backward 
differences taken at n times the step of A and V. A relation between and A is 
then given by 


(l-V n ) = <l+A)-» 
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Hill II I II 


since 


(l- v n ) f ( x ) ~ f < x “ **) 


= E" n f(x) 

= (1 + A)~ n f (x) . 


Analogous to (4) , 


A = (1-V )-V» - ! , 


(IS) 


and the multirevolution predictor or extrapolation formula can be derived 
by expressing the forward difference A n in terms of backward differences vj of 
the forward difference A . From (1) and (15) , in analogy with (5) , 


A 

n 


v 

n 




V 


n 


l 

(!— V n)' Vn -l 


A 



(16) 


where the coefficients y , which are polynomials in 1/n, can be determined using 
the series (6) and the expansion 

co 

(1 " X)a = IW (i) xi • (17) 

i= 0 

We have 

00 

(i-v„r”-i = £i-iy CY n ) v „‘ 

i=i 
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IIIMHII IIIIINMI Hill llll I 


where b Q = 1 and 




b k = (-l) k 


j=l 


(k+ryr 


k = 1,2, 3, 


(18) 


Hence 


00 

" !>» 


i = 0 


1 - V, 


- (1 - v ,')' I/n - 1 


oo 

E V 

i=0 


00 

= n y i V n i ’ 


where the y i 's are given recursively . y 


i = 0 


7n “ 1 


and 


y, 


E b J n-. 


j=i 


Using (16) and (19), we then obtain 

A n f ( X ) = f ( X +nh ) “ f ( x ) 


= + (i-s) v n + n(s-| 




r < 19 ) 


1, 2, 3, ••• 




+ • ,, |nAf(x) , (20) 
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which is the formula used to compute the predicted or extrapolated value of 
f(x +nh) by using the backward differences V n ‘ of the forward difference of f(x) 
at l/n th of the step. 

Again, in analogy with the integration formulas, the corresponding multi - 
revolution corrector can be derived by expressing the backward difference in 
terms of backward differences of A. 


The result is (see (11)) 


V 

n 


v 

n 




A 



( 21 ) 


where, as above, we find that the coefficients y.* are given recursively by 

1 


7n* 


~\ 


and 




* = 


1 

- £ b i 7 *i 


> ( 22 ) 


i = 1, 2, 3, • 


j=i 




where the bj 's are given by (18). 
Using (21) and (22), 


V n f ( x ) - f ( x ) “ f(x -nh) 

= (l + Kk-T7(l-^H 2 - 


•••}"* f ( x ) • (23) 


which is the formula used to compute a corrected value of f(x) by using a for- 
ward difference of an extrapolated value for f(x). 



The similarity between the coefficients given by (8) and (19) and (12) 
and (22) should be noted, as should, in analogy with (13) , the relation 





Finally, in practice, an algebraic equivalent of Equations (10) and (13), 
known as the "summed" form of the difference formulas, is often used. It has 
been established (Henrici, 1962) that this modification considerably reduces 
the propagated roundoff error associated with these formulas. Formally, one 
can obtain this modification by applying the inverse operator V ' 1 to both sides 
of Equations (10) and (13). In like manner, one can obtain the summed form of 
the extrapolation formulas by applying the operator V n _1 to both sides of 
Equations (20) and (23). 


THE MULTIREVOLUTION ALGORITHM 

Let f . denote the value of an orbital element at the descending node of the 
j th revolution , and let n be the number of revolutions to be stepped. Also let 
k be the order of the highest backward difference V 1 to be retained in formulas 
(16) and (21). We then have the multirevolution formulas 


and 



(24) 


(25) 


We see that, in analogy with the Adams integration process, a set of start- 
ing values must be computed before these formulas can be used. Assume then 
that the values Af . = f , + 1 - f . , i = 0, n, 2n, 3n, • • • , kn have been computed, 
i.e., we assume knowledge of the orbital elements at descending nodes 0, 1, n, 
n + 1, 2n, 2n + 1, • • • kn, kn + 1. Using these values, we can compute the 
differences 

V n ( Af k n) i = 1, 2, ••• k , 
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The multirevolution cycle consists of the following steps: 

(1) Extrapolate the orbital elements n revolutions ahead using (24) with 
j = nk, i.e., compute a predicted value for f (k+ 1)n , the element at the de- 
scending node of the (k + l)n th revolution. 

(2) Starting with these extrapolated values, integrate one revolution to 
obtain the elements at descending node (k + l)n + 1 and form the forward difference 

^^(k + l)n ~ f(k+l)n+l ' f(k+l)n ' 


t 




(3) Form the backward differences^ 1 (Af (k+1)n ) , i = 1, 2, . . . k, and use 
Formula (25) with j = (k + 1 )n to obtain a corrected value for f (k+1)n • The 
cycle is complete. Repeat the process with j = (k + l)n in step (1) . 

The following remarks can be made concerning the computation procedure: 

(1) Values of the orbital elements at the descending node of any revolution 
can be obtained by inverse interpolation, i.e., first, a direct interpolation is 
made to find the time of nodal crossing, followed by an inverse interpolation to 
obtain the values of the orbital elements. 

(2) Other points in the orbit such as perigee, apogee, or ascending node 
could be used instead of the descending node. 

(3) The time associated with the extrapolated orbital elements can also be 
obtained by the extrapolation process, allowing one to obtain values for the orbital 
elements at any point in a computed revolution. 

(4) For any particular orbit, the number of revolutions that can be 
stepped by this process and yet maintain some predictable accuracy can be 
determined by comparing the change in the orbital elements from node to node 
with the change required by the short-step integration. 

(5) For any particular orbit, the number of differences V n ‘ to be retained 
in the extrapolation formulas can be determined precisely as in the case of the 
integration formula, i.e., retaining differences of sufficiently high order so that 
the "truncated" terms are within the accuracy of the computations. 


NUMERICAL RESULTS 

The equations of motion considered in the numerical computations are of the 
form 

- /xx — — 

x = - ^3 + F 1 + F 2 > (26) 
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where x = (x, y, z), R = (x 2 +y 2 + z 2 ) 1/2 , Fj and F 2 are perturbation functions, 
and fi is a constant. The formulation of these perturbations and the values of 
the related constants used are given in the appendix. 

Before beginning, the following remarks concerning the numerical 
results can be made: 

(1) All short -step numerical integrations were performed with an Adams - 
Cowell 13th-order predictor-corrector process (Henrici, pp. 192-194, 291- 
293), using a Runge-Kutta-type starter, so that this process was used both for 
the comparison integrations (Table 1) and for the short -step integrations 
required by the multirevolution algorithm . 

(2) A predictor -corrector tolerance of 1 x 10' 12 was used for all the 
Adams -Cowell integrations. The number of derivative evaluations tabulated 
was obtained by multiplying the number of steps taken at the Cowell time step h 
by the average number of predictor -corrector iterations taken over the range 
of integration. 

(3) For the examples, the integrations were performed for approximately 
100 revolutions of three orbits. 


Initial Elements* 

Example I 

Example II 

Example III 

Semimajor axis (Earth radii) 

1.15 

1.26 

6.71 

Eccentricity 

.075 

.072 

.003 

Inclination 

1.52 

1.03 

.0004 

Mean anomaly 

6.03 

3.71 

3.80 

Argument of perigee 

1.15 

3.14 

.31 

Long, of ascending node 

4.76 

6.16 

2.29 

Period (min) 

104.25 

120.37 

1469.46 

Range of integration (min) 

11,000 

12,500 

148,000 


(4) All computations were performed on the UNIVAC 1108 computer using 
double-precision arithmetic. 

In Table 1, the results of the comparison integrations are tabulated. Error 
estimates were obtained by integrating Equations (26) with F t = F 2 = 0 and 
comparing with the Keplerian solution. For the cases considered, the approximate 
ratio of the change in the position that occurred at the Cowell time step h and 
the corresponding change that occurred from descending node to descending 
node were found to be 10 or greater. 


•Angles in radians. 
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Table 1 


Orbit 

h 

(min) 

No. of Rev. 
Computed 

Total No. of 
Der. Eval. 

Estimated 

Error 

Computation Time 
(min) 

I 

1 . 0 

105 

17,337 

2 x 10~ 9 

1.14 

II 

.8 

104 

15,617 

4 x 10" 12 

1. 09 

III 

22. 0 

100 

7,257 

9 x 10~ 9 

. 71 


Consider the results obtained by applying the multirevolution algorithm to our 
test cases. The elements used in the extrapolation process were the position and 
velocity vectors themselves. (The time associated with the extrapolated values was 
also obtained in this manner. ) We begin by tabulating the results obtained by using 
predictor formula (24) only (Table 2). The error estimate given is the difference 
between the position at the node of the 100th revolution obtained by the comparison 
integration (Table 1) and by the multi revolution process. In each case, for each n, 
the order of the process (k) was increased until no change in accuracy occurred. 


Table 2 


Orbit 

n 

k 

Total No. of 
Der. Eval. 

Estimated 

Error 

Computation Time 
(min) 

I 

3 

2 

6049 

1 x 10' 8 

0.55 



4 

6737 

3 x 10' 9 

.59 


5 

2 

4586 

3 x 10~ 8 

.39 



4 

5934 

2 x 10' 9 

.48 


7 

2 

4337 

8 x 10" 8 

.35 



4 

6343 

4 x 10" 9 

.48 


9 

2 

4542 

1 x 10' 7 

.36 



4 

7209 

1 x 10' 9 

.52 

II 

3 

2 

5650 

9 x 10' 8 

.51 



4 

6268 

8 x 10' 11 

.55 


5 

2 

4263 

5 x 10' 7 

.36 



4 

5483 

6 x 10~ n 

.44 


7 

2 

4013 

1 x 10' 6 

.33 



4 

5834 

4 x 10" 10 

.45 



6 

7656 

8 x 10" n 

.57 


9 

2 

4188 

3 x 10" 6 

.33 



4 

6611 

2 x 10' 9 

.49 



6 

9034 

4 x 10~ 10 

.65 

III 

3 

2 

2329 

6 x 10' 7 

.62 


5 

2 

1803 

5 x 10" 7 

.45 


7 

2 

1815 

5 x 10' 7 

.40 


9 

2 

1923 

5 x 10“ 7 

.37 
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An examination of these results shows that significant gains in efficiency were 
obtained by application of the algorithm to the numerical integrations, with no 
significant loss in accuracy. 

A good portion of the computation time tabulated was used in forming the 
starting table (Figure 1) and in restarting the short-step integrations required 
by the process, since these were done with the Runge-Kutta, Adams-Cowell 
process. 

In Table 3, the results obtained by applying the predictor formula (24), 
followed by an application of the corrector formula (25) , are tabulated. As in 
Table 2, for each n , the order (k) was increased until no change in accuracy 
occurred. 


Table 3 


Orbit 

n 

k 

Total No. of 

Estimated 

Computation Time 


Der. Eval. 

Error 

(min) 

I 

3 

2 

10,903 

2 x 10” 9 

0.99 


5 

2 

7,317 

5 x 10' 9 

.65 



4 

8,361 

3 x 10' 9 

.70 


7 

2 

6,157 

1 x 10" 8 

.52 



4 

7,860 

1 x 10 -9 

.62 


9 

2 

5,907 

2 x 10~ 8 

.48 



4 

8,271 

2 x 10” 9 

.62 

II 

3 

2 

10,194 

2 x 10” 8 

.92 



4 

10,528 

8 x 10” u 

.94 


5 

2 

6,819 

7 x 10" 8 

.59 



4 

7,755 

7 x 10” 11 

.65 


7 

2 

5,717 

2 x 10 -7 

.48 



4 

7,254 

6 x 10' 11 

.58 


9 

2 

5,466 

4 x 10‘ 7 

.44 



4 

7,605 

1 x 10 _1 ° 

.58 



6 

9,744 

2 x 10 -11 

.72 

in 

3 

2 

4,086 

6 x 10" 7 

.99 


5 

2 

2,781 

5 x 10” 7 

.68 


7 

2 

2,442 

5 x 10" 7 

.56 


9 

2 

2,387 

5 x 10 -7 

.48 


These results show that an application of the corrector yielded slight gains 
in accuracy at some expense in efficiency. Again, as in Table 2, this cost in 
efficiency was due to the process used for all the short-step integration required 
by the algorithm. 
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At this point, several suggestions should be considered. 

(1) As suggested by Cohen and Hubbard and Mace and Thomas, an iterative 
process using the multirevolution formulas, instead of the short-step integration 
process, could be used to form the required starting table (Figure 1). Such a 
scheme would increase the efficiency of both the predictor and predictor -corrector 
processes used to obtain the results. 

(2) Since much of the computation time, especially in Table 3, was used in 
restarting the short-step integrations by the Runge-Kutta technique, an increase 
in efficiency could be obtained by extrapolating the required values for these 
integrations, as well as the nodal points. 

(3) With more efficient starting processes, it could be feasible to iterate on 
the corrector formula (25) to improve the accuracy when using large n . This 
would be analogous to the Adams integration process. 


SUMMARY AND CONCLUSIONS 

The multirevolution algorithm has been reviewed and found to be an effective 
tool applicable to the orbit determination problem. 

The recursive formulas obtained for the coefficients required by the process 
have been found to be an effective and simple means of obtaining these coefficients 
for any given n and k. The numerical results have indicated that the "predictor 
only" scheme can be applied to the common multistep integration process with 
a Runge-Kutta starter to yield considerable gains in efficiency and, if more ef- 
ficient starting processes were used, gains in accuracy due to an application of 
the corrector could be obtained at little or no expense in efficiency. 


Goddard Space Flight Center 
National Aeronautics and Space Administration 
Greenbelt, Maryland, January 29, 1970 
311-07-21-0 1-51 
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Appendix 

THE EQUATIONS OF MOTION 


The equations of motion considered in this paper are of the form 


and 


ux — — 

x = - iT + F 1 + F 2 ’ 


(to) 


*(t 0 ) 


X 0 ’ 


where x - (x, y, z ), R = I x | = ( x 2 + y 2 + z 2 ) 1/2 , ,u_is a constant, x 0 and i? 0 
are the given initial position and velocity vectors, Fj is the perturbation due 
to the nonsphericity of the Earth, and F 2 is the perturbation due to drag. 


F 


1 



is given by 

F x = § [jR 2 (5s-l) + Hz(7s-3) + | (-63s 2 + 42s - 3 )] , 

F y = ^ [jR 2 (5s-l) + Hz(7s-3) + f (-63s 2 + 42s - 3 )] , 
and 

F z = [j R 2 (5s - 1) + ^- (- 63s 2 + 70s - 15)] + (35s 2 - 30s + 3 ) » 
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where s = (z/R) 2 and J,H, and K are the second, third, and fourth harmonics 
of the Earth's potential field. F 2 * s °f the form 


C D A , _ . _ 

2M V r | V r ’ 


where V r is the relative velocity of the satellite, p the atmospheric density* at 
the satellite position, andA,M, and C D are the cross-sectional area, mass, and 
drag coefficient of the satellite, respectively. 

Numerical constants: 

Unit of Time = 806.832 sec 
Unit of Length = 6378.388 km 

p. = 1 

J = +0.162 x 10~ 2 
H = -0.640 x 10" 5 
K = +0.690 x 10" 5 
A = +0.116 x 10 5 cm 2 
M = +0.158 x 10 6 gm 
C D = +2.3 


•The density was computed (from a standard model) for satellite heights of less than 1000 km; otherwise p 
was set to zero. Hence, in our examples, only orbit I was influenced by F 2 * 
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